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ANALYTICAL  HAZARD  REPRESENTATIONS  FOR  USE 
IN  RELIABILITY,  MORTALITY,  AND  SIMULATION  STUDIES 

Part  I 
D.  P.  Gaver  * 
M.  Acar 


I.   INTRODUCTION 

The  failure  rate  function,  or  hazard  function  (hazard 
for  short)  may  be  described  as  the  conditional  probability  of  an 
equipment's  failing  at  operating  age   t,  having  survived  to  that 
age.  The  reliabilities  of  a  variety  of  electronic  and  mechanical 
items  are  conveniently  and  naturally  described  in  terms  of  the 
appropriate  hazard  function,  and  so  is  the  longevity  of  human 
beings.   The  term  force  of  mortality  replaces  hazard  in  the 
latter  context. 

This  paper  is  devoted  to  a  study  of  several  simple 
analytical  representations  for  hazard  functions.   These  repre- 
sentations are  in  turn  based  upon  representations  of  random 
variables  having  certain  required  properties,  in  terms  of  others 
having  familiar  distributions — in  particular  the  exponential. 
Similar  ideas  are  due  to  Tukey  (19  76)  and  recently  have  been  examined 
by  Parzen  (1973).  The  hazard  representations  proposed  are  quite 
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expeditiously  used  in  simulation  studies,  e.g.  of  system 
reliability  or  availability  in  terms  of  component  lifetimes. 
They  may  also  be  used  in  data  analysis  studies,  in  order  to 
parsimoniously  describe  data  sets  in  terms  of  perturbations  of 
convenient  and  familiar  standard  distributions.   Their  use  in 
data  analysis  and  simulation  is  also  described  in  Gaver,  Lavenberg, 
and  Price  (1976),  and  in  Gaver  and  Chu  (1977). 


II.   SYSTEM  FAILURE  PATTERNS 

It  is  plausible  to  think  that  the  time  series  of  failures 
in  a  system  may  involve  these  stages. 

1.  Early  failures.   There  may  be  a  relatively  large  number  of 
failure  soon  after  a  system  is  introduced  because  of  design 
defects,  production  errors,  or  errors  stemming  from 
maintenance  personnel  inexperience.   This  situation  is 
characterized  by  a  hazard  function  that  is  initially  large, 
but  that  decreases  with  time.   "Infant  mortality"  is  in 
evidence . 

2.  Random  Failures.   Following  the  early  failure  period  there 
may  be  a  period  during  which  failures  occur  at  an  essentially 
constant  rate  for  a  rather  prolonged  time.   During  this 
period  the  hazard  function  is  nearly  constant,  so  the  times 
between  failures  are  close  to  being  exponentially  distributed 
The  effect  of  age  or  wearout  is  not  yet  apparent. 


3.   Wearout  Failures.   Eventually  following  the  period  during 

which  a  constant  hazard  is  evident  there  is  likely  to  be 

a  period  of  ever-increasing  failure  rate  caused  by  wearout 
of  system  components. 

A  graphical  representation  of  a  hazard  function  that 
exhibits  the  behavior  described  is  given  below.   Note  that  it 
has  the  legendary  "bath  tub"  shape. 


Hazard 


time 


Some  comments  on  the  above  follow: 

A.  The  term  "failure"  may  refer  to  an  event  that  is  analogous 
to  human  death,  after  which  the  entire  system  is  replaced. 
On  the  other  hand  repair  or  component  replacement  may  occur 
after  failure;  the  system  is  only  repaired,  not  entirely 
replaced.   In  the  former  case,  a  hazard  function  of  the 
kind  depicted  in  Figure  1  applies  to  each  system  event 

("death");  when  the  system  is  installed  (or  is  born)  that 

hazard  operates  starting  from  scratch  at   t  =  0   until 

system  failure  (human  death,  for  instance)  after  which  a 

similar  hazard  goes  into  effect,  starting  once  again  from 

zero.   In  the  latter  case,  in  which  repair  of  a  component 

occurs,  a  hazard  function  like  that  of  Figure  1  applies  at 

t  =  0,  but  after  the  first  event  ("failure")  at   t, ,  a 

repair  action  is  accomplished.   The  same  hazard  operates 

for   t  >_  t.    until  the  next  event  at   t?  >  t,  ,  and  so 

on.   Intermediate  situations  may  be  envisioned,  in  which 

after  event  n   at   t   the  hazard  governing  system  failure 

n+1   starts  at   t   -t,0<t   <t. 

n    n       n    n 

B.  Although  there  is  reason  to  assume  that  hazards  somewhat 
like  that  of  Figure  1  occur  in  general  for  systems,  the 
possibility  exists  that  the  system  hazard  is  "bumpy" 
because  wearout  failures  of  components  or  subsystems 
may  well  occur  at  intermediate  times. 


C.   If  the  theory  is  applied  to  systems  with  little  or  no 
wearout  propensity,  as  should  be  the  case  when  dealing 
with  computer  software  modules,  then  the  hazard  function 
may  well  exhibit  the  initial  falloff  of  Figure  1  but  not 
the  rise  at  later  times.   In  fact,  a  constant  decline  as 
bugs  are  found  and  removed  could  be  (optimistically) 
anticipated  for  software.   The  right-hand  side  of  the  bath 
tub  vanishes,  and  the  picture  is  that  of  a  ski  slope. 


III.   MODELS  FOR  THE  HAZARD  FUNCTION 

In  this  section  mathematical  models  are  presented  for  the 
failure  rate  or  hazard  function.   Recall  that  the  hazard  may  be 
defined  as  follows. 

Definition.   Suppose  that  the  time  to  failure,  X,  is  a  random 
variable  with  distribution  function   F(x),  where   F(0)  =  0; 
the  latter  possesses  the  density  function   f(x),  f(x)  =  dF/dx , 
such  that  for  any  positive   x, 


x 
F(x)  =  /   f(y)dy  .  (3.1) 

0 


Then  the  hazard  function,  or  failure  rate  at  age   x,  is  given  by 


h(*»  "!-=%.  <3-2' 


The  interpretation  of   h(x)dx   is  that  it  is  the  conditional 
probability  of  failure  in  the  interval   (x,  x  +  dx) ,  given 
that  there  has  been  no  failure  up  to  age   x. 

Express  the  hazard  as 


,    ,    N  dF/dx 

h(x)    = 


1    -    F(x) 


h(x)dx   =    1    _dp(x)    =    "   d{log[l    -   F(x)]}, 


it  then  follows  after  integration  that 


x 
F(x)  =  1  -  exp[-  /   h(y)dy]  .  (3.3) 

0 


Thus  if  the  hazard  is  specified,  so  is  the  distribution  function, 
and  conversely. 
Note  that  if 

h(x)=A>0,  0<_x 

then 

F(x)  =  1  -  e"Ax,  0  <  x  ,  (3.4) 

so  a  constant  hazard  function  implies  the  exponential  distribution 
of  the  random  variable   X,  and  conversely. 


Obviously  a  constant  hazard  representation  does  not 
describe  the  bath  tub  hazard  shape  of  Figure  1,  nor  does  it 
represent  a  situation  in  which  hazards  decline,  possibly  because 
design  defects  or  "bugs"  are  occasionally  removed.   Here  are  two 
hazard  representations  likely  to  be  useful  for  such  purposes. 

(1)   A  Bath  tub  Model 

Define  the  random  variable   Z   in  terms  of   X,  X   being 
exponentially  distributed  with  mean   A   ,  as  follows: 

Z  =G(X)  =  XL(X)  R(X) 
or  (3.5) 

=  XcJ)  (X)  ,     <J>  (X)  =  L(X)R(X) 

where 

a)  L(x)   is  convex  in   x,  L(0)  <  1,  L(°°)  =  1, 

b)  R(x)   is  concave  in  x,  R(0)  =  1,  R(0)  >  R(°°). 

Then  the  hazard  of   Z   may  be  made  to  exhibit  a  bath  tub  shape, 
as  in  Figure  1,  by  proper  choice  of  the  function   L   and   R. 

Example.   Suppose 

L(x)  =  i  +*  v  a  >  0,  0  <  x 

1  +  ax 

(3.6) 


R(x)  ■  rrV         e  >  °  ' 


Clearly 


■r  /  \  T-.  /  \        Otx         1        ax         x 
z  =  xL(x)  R(x)  =  x 


1  +  ax   1  +  3x   1  +  ax   1  +  3x 

is  a  monotonically  increasing  function  of   x.   Furthermore, 

_3 

choose   a   large  (e.g.  a  =  10)  and   $   small  (e.g.   3  =  10   ). 

Then  it  is  intuitively  clear  that   (i)  small  x-values  transform 
into  even  smaller  z-values,  e.g.   x  =  1   corresponds  to   z  =  0.91 
and  x  =  2   corresponds  to   z  =  1.90,  but   (ii)  this  effect 
dwindles  as   x   increases,  so   x  =  10   corresponds  to   z  =  9.8 
and   x  =  50  to   z  =  47.5   and  the  z-values  closely  resemble 

the   x's  percentage-wise,  but   (iii)  as  x   increases  still 

3 

further  the   z's  do  not  follow  suit:   x  =  10    corresponds  to 

z  =  500.   This  suggests  that  if   x   is  a  value  assumed  by   X, 
that   Z   shares  the  properties  of   X   in  mid-range,  i.e.  for 
intermediate  x-values,  but  differs  from  X   by  having  a  dispro- 
portionate probability  of  assuming  small  values  (near  zero) , 
or  large  values  (near,  but  less  than,  1/3) .   Thus  the  hazard  of 
Z   will  appear  to  be  a  "bath  tubbed"  version  of   X,  particularly 
if   X   is  exponential. 

We  focus  attention  on  the  representation  (3.6)  in  what 
follows,  mainly  for  analytical  and  computational  convenience. 
Of  course  there  are  many  other  possibilities,  such  as 


L(x)  =  1  -  e 


R(x)  =  e3x 


(3.7) 


these  latter  may  be  adjusted  to  provide  sharper-edged  bath  tubs 
than  can  (3.6),  but  iteration  of  (3.6)  may  be  induced  to 
accomplish  the  same  purpose. 

( 2 )   A  Decreasing  Failure  Rate  Model 

Define  the  random  variable   W   in  terms  of   X,  X   again 
being  exponential  with  parameter   A   : 

W  =  XT(X)  (3.8) 

where   T(x)   is  an  increasing  function  of   x,  L(0)  =  1.   Then 
the  hazard  may  be  made  to  exhibit  a  decreasing  behavior. 

Example .   Suppose 

T(x)  =  1  +  ex  ,         c  >  0,  0  £  x.  (3.9) 

Then 

z  =  x(l  +  ex)  (3.10) 

is  monotonic,  and  small  x-values  lead  to  comparable  z-values 
(especially  when   c   is  small) .  but  larger  x-values  are  "amplified" 
by   1  +  ex   to  yield  increasingly  large  z-values. 


Attention  will  be  focused  upon  (3.9),  although  other 
possibilities  exist  that  accomplish  the  same  purpose,  namely 
that  of  lengthening  the  right  tail  of  the  distribution  of   X 
(simulating  outliers,  for  instance)  while  leaving  the  body  of 
the  distribution  virtually  unchanged. 


IV.   MATHEMATICAL  PROPERTIES  OF  THE  "BATH  TUB"  HAZARD  MODEL 

Various  analytical  properties  of  the  previously  described 
models  will  now  be  recorded.   These  provide  useful  insights  into 
the  behavior  of  the  random  variables   Z   and  the  underlying 
(generating)  variables   X. 

A.   Monotonicity;  Quantiles 

It  is  convenient  to  focus  on  montonic  increasing  trans- 
formations, i.e.  if 

z  =  G(z)  =  z<J>(z)  (4.1) 

then  in  order  that  the  above  function  be  monotonically  increasing, 

dz/dx  >  0.   Observe  that  logarithmic  differentiation  of  (4.1) 

provides 

dz=dx+1|(x) 
z    x     <j>(x) 

and  thus   dz/dx  >  0   if  and  only  if 
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1  +  *1M     >  o  (4.3) 

x     d>  (x) 


Alternatively,  the  condition  is,  in  terms  of   L(x)   and   R(x), 

1    LMK)     RMx) 
x     L(x)     R(x) 


It  is  easily  seen  that  the  important  example  (3.6) , 


i    \  ctx 

(x)  = 


1  +  ax    1  +  $x  ' 

yields  a  monotonic  relationship  between   z   and  x.   The  fact 
that  this  transformation  can  be  easily  and  explicitly  inverted 
(solved  for   x   in  terms  of   z)  will  be  exploited  subsequently. 
Of  course  if   z(x)   is  monotonically  increasing  then  so 
is   x(z),  the  inverse  function.   The  events   (Z  <  z)   and 
(X  _<  x(z))   are  equivalent,  and  so 

P{Z  <_  z}  =  P{X  <_  x(z)  },  (4.5) 

from  which  it  follows  that  if   x   e  x (p)   is  the  p-100%  quantile 
of   X,  i.e. 

P{X  <  x(p) >  =  p  ,  (4.6) 

then 

P{Z  <  z(p)}  =  P{Z  <  z(x(p))}  =  p  (4.7) 
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and  so   z(p),  the  p-100%  quantile  of   Z   is  simply  obtained  from 

z(p)  =  x(p)  <Mx(p))  =  x(p)  L(x(p))  R(x(p))      (4.8) 

In  other  words  we  very  easily  translate  from  (points  on)  the 
inverse  distribution  of   X   to  the  inverse  distribution  of   Z. 
Explicit  representation  of  the  distribution  of   Z   is  however, 
not  often  easily  possible. 


B.   Hazard  and  Density  Function  Relationships 

In  order  to  investigate  the  relationship  between  the 
hazards  of   Z   and   X,  begin  by  writing 

x(p) 

0 


p  =  Fx(x(p))  =  1  -  exp[-  /     h  (u)du]  (4.9) 


or 

x(p) 
/     h  (u)du  =  -  £n(l-p) 
0      x 


Now  differentiate  with  respect  to   p   to  find 


x        dp 
or 


<x<p»^  =  r^  (4.io, 


hx(x(P,)  '&'   iVV*'^  •  iV       (4-11> 

nere   hx   and   f^   are  the  hazard  and  density  functions  of  the 
r.v.   X.   The  relationship  (4.11)  holds  for  any  distribution,  of 
course . 
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Differentiation  of  (3.5)  reveals  the  connection  between 

h    and  h  .   From  (4.2) 
z        x 


dz(p)     /  , 
-3p   =  z(p) 


1   ,   <j>'  (x(p)) 
xTpT    <j>(x(p)) 


dx(P; 
dp 


(4.12) 


=  U(x(p))  +  x(p)  <j>'(x(p))]  ^|2l 
From  (4.11),  applied  now  to  the  z-hazard,  there  results 


h  (z(p))  =  h  (x(p))  [*(X(P»)  +  X(P>  *'(x(p))]  ,  (4.13) 

Z  X 


so 


h  (z(p))  =  h  (x(p))  .,    i    t  *    _,_ ,\    tt-7  /  i  >      (4.14) 

zv  ^      x   ^   4>(x(p))  +  x(p)  4>'(x(p)) 

Multiplication  of  both  sides  by   1-p   then  shows,  in  view  of 
(4.11),  that  the  density  functions  are  similarly  related: 

f2<z(p>)  =  fx(x(p))  ^x(p))  +  K)p)    „.(x(p>) 


Example. 

X   is  exponential (X) .   Then 


hz(V  "  <|>(x(p))  +  x(p)  <j>'  (x(p)) 


(4.15) 
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Now   use    the   specific      <p  (x)      of    (3.6): 


4>(x)     =    j-+ 


ax        1    +    3x 
or,    in   terms    of   logarithms, 

Zn    4>(x)    =    In   ax   -    £n(l    +   ax)    -    £n(l    +    3x)     , 


so 


4>'  (x)         1  a  3_       .  1   -   a3x 

cf)(x)         x        1   +   ax        1   +   3x        x(l   +   ax)  (1   +    3x)     ' 


(4.16) 


and 


<J>  (x)    +   X<J>'  (x)    =    <|>(x) 


2   +    (a   +    3)x 
(1   +   ax)  (1   +    0x)" 


(4.17) 


finally 


h    fzfnn    -   Ml   +   ax(p))2    (1   +   3x(p)) 
nzlziP"  ax(p)     [2    +    (a   +    3)    x  ( p7 


PIT 


(4.18) 


Although  this  expression  is  not  quite  explicit,  qualitative  properties 


of   h    can  be  deduced  from  it. 
z 


(a)   if   p  -*>  0 ,  x(p)  =  -  y  £n(l-p)  i    0,  and  hence 


*,<*<P>>~  2H"W  ' 


(4.19) 


or 


lim  x(p)  h  (z(p))  =  £ 

p  -»■  0 


(4.20) 
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Since  for   p  -*•  0 , 


and  hence 


there  results 


or 


2 

z  (p)  ~  ax  (p) 


x(p)  ~  [z(p)/a]1/2  (4.21) 


h   (Z(P))   ~  ,-rj 

Z  2[z(p)]1/^  /a" 


lim    ZzTpT  h  (z(p))  =  — - —  (4.22) 

p  ->  0         z         2  /a 


This  shows  that   h  (z)  +  °°   as   z  ■*■  0 ,  creating  the  left-hand  end 


of  the  bath  tub  of  Figure  1. 


(b)   If   p  -*■  1,  x(p)  \    °°,  and   z  (p)  i    1/8   so 


2   2 


„,(,<*»-  xHiT_x_(|) 


or 


For   p  ■>  1 


so 


and  thus 


lim    r777?  V2(p>>  =  x  ^  (4-23> 

p    ->    °°     [x(p)  ] 


1    -    8z(d)  ~  ilL-t-J-i (4    24) 

X         PZip;       a3      x(p)  K*-**) 


*<p>  ~  ^  i^w  (4-24) 


15 


a  +  e 

a 


h„(z(p))    X\        i        J  — j  (4.26) 


(1  -  3z) 


Once  again  it  appears  that  the  hazard  rises  rapidly,  this  time 
as   x(p)  +  °°   and   z  (p)  i    3   ;  the  other  end  of  the  bath  tub 
is  thus  fashioned. 


(c)   If   p  =  1  -  e~  ,  then   x(p)  =  A~  .   Then 


h^zd-e"1))  -  Ul  +  a/A]2  [1  +  3/A]2  (4>2?) 

j      [2  +  (a  +  3)/A] 


The  bath  tub  effect  is  presumably  achieved  by  choosing   a   large 
and   3   small.   Let   a  •*■  0   and    3^0   independently  in  (4.27)  ; 
it  is  clear  that  the  limiting  value  of  the  hazard  is   A.   This 
indicates  that  the  hazard  is  (approximately)   A   for  middling 
values  of   z. 

C.   An  Explicit  Formula  for  a  Hazard 

The  expression  (3.6)  leads  to  the  relationship 

???,(P),  ...^  ,  (4.28) 


Z(P)  "  (1  +  ax(p))  (1  +  $x(p)) 


and  the  latter  may  be  explicitly  inverted  by  solving  a   quadratic 
equation.   The  result  is 
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x(p) 


(a  +  3)  z(p)  +.  /(a  +  3)2  z2(p)  +  4z(p)  a(l-3z(p)) 

2a (1  -  3z(p)) 


(4.29) 


Wow  a  direct  differentiation  of  this  expression  and  invocation 
of  (4.11)  produces  the  expression 


hz(z)  "  2a(l-3z) 


3((a+3)z  +  /(a+3)2z2  +  4az(l-3z)  }" 


1  -  3z 


+  (a  +  3) 


{(a-3)2z    +   2g}    /(g+3)2z2    +   4~az  (l-3z) 
(a   +    3)2z2    +   4az(l    -    3z) 


(4.30) 


This  form,  while  explicit,  provides  no  particularly  useful  insights; 
the  bath  tub  end  shapes  already  noted  in  (4.22)  and  (4.26)  can  be 
deduced  directly  from  (4.30). 

Some  graphical  plots  of   h    are  presented  later.   They 
illustrate  the  behavior  of  the  present  hazard  representation  in 
a  more  understandable  fashion  than  does  the  formula  itself. 

D.   An  Explicit  Formula  for  the  Failure  Time  Distribution 

Because   x   and   z   are  monotonically  related  through 
(4.28)  we  have 
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P{Z    <    z}    =    P 


)   Y    <    v    _     (a+3)2    +    /(g+e)^z2    +    4gz(l-3z)       I 
|   X   -  X  2a  (1    -    32) 

=   1    -   exp  )  -A       (a+3)z    +    /(a+3)2z2    +   4az(l-3~z7     1  J 


exp  /  -A 


2a (1  -  3z) 


(4.31) 


Again  the  explicit  formula  seems  unproductive  of  insights 


V.   MATHEMATICAL  PROPERTIES  OF  THE  DECREASING  FAILURE  RATE  MODEL 
A.   The  Hazard  Behavior. 

The  expression  (4.14)  can  be  applied  to  deduce  the  hazard 
function  of  the  representation  (3.8),  advertised  to  produce  a 
decreasing  failure  rate.   There  we  specified 


<J)(x)  =  T(x)  =  1  +  ex  , 


(5.1) 


and  thus,  from  (4.14), 


hw(w(p))    (1  +  cx(p))  +  x(p)c  "  1  +  2cx(p) 


(5.2) 


Qualitative  properties  follow  easily. 


(a)   If   p  ■>  0,  x(p)  \    0,  and 


h  (w  (p)  )  ~  A 
w   ^ 


(5.3) 


Thus  the  hazard  is  approximately   A   for  small   z 
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(b)   if   p  ■*■   1,  x(p)  t  °°,  and 

2 
w(p)  ~  c[x(p)  ] 

so 

h  (w(p))  ~      A      ,  (5.5) 

2/cwTpT 

which  clearly  decreases,  as  claimed.   It  may  be  inferred  that  the 

distribution  of    W   appears  nearly  exponential,  but  possesses 

an  extraordinarily  long  right  tail--possibly  the  result  of  outliers 


B.   Explicit  Formulas  for  the  Hazard  and  the  Distribution  Function 
Direct  solution  of  the  quadratic  equation 


2 
w  =  x(l  +  ex)  =  x  +  ex 


presents 


x(p,  =  "  +  4^(p)  -  1  f 


(5.6) 


which,  when  differentiated,  leads  to 


h    (w)    = h    (x) 


w  /1~T4 


cw 


X 
/l    +    4cw 


(5.7) 
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The  distribution  function  is 


v  s  /l   +   4CW  -   1  (  ,r       n\ 

P{W<w}    =   P{X<_x  =   * j  (5.8) 


=    1    -   exp  <  -  A 


/l   +    4cw   -    1 
2~c" 


This  distribution  bears  a  close  family  resemblance  to  the  Weibull 
distribution   1  -  F(w)  =  exp{-k/w},  especially  for  large  (right 
tail)  values  of  w. 


VI.   AN  ALTERNATIVE  "BATHTUB"  HAZARD  REPRESENTATION 

The  simple  parametric  model  (3.6)  leading  to  a  bathtub-shaped 
hazard  is  by  no  means  the  only  possibility.   We  next  describe 
another  simple  approach.   It  is  that  of  defining  a  hazard  function 
having  an  appropriate  shape ,  and  then  deducing  the  corresponding 
distribution  function,  and  a  procedure  for  sampling  from  it, 
rather  than  proceeding  in  reverse  order,  as  before. 

Let  the  hazard  be  of  the  form 

h(z)  =  g(z)  +  A  +  k(z)  ,  (6.1) 


where   g(z)  >  0   is  a  decreasing  function  of   z   such  that 
lim  ^  ^  g(z)  =  0,  and  k(x)   is  an  increasing  function  of   x, 
such  that  (preferably)  k(0)  =0   and  k  (°°)  =  °°.   Such  a  function 
can  yield  a  bathtubbed  hazard. 
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Example 


h(z)  =  z  |  g  +  Bz  +  X  (6.2) 


A,  B,  a,  A   all  positive. 

Clearly,  (6.2)  has  a  generally  "bath  tub-like"  appearance, 
since 


for  if 


h'  (z)  =  -  - 7   +  B  (6.3) 

(z  +  a) 


-  -  +  B  <  0,      then   h' (0)  <  0  (6.4) 


while  for 


■d)1 


z  >  zn   =  I  *)  _  a  (6.5) 


h(z)  >  0. 

Detailed  behavior  is  adjustable  by  choice  of  the  parameters 
Now  the  distribution  function  of  time  to  failure,  Z,  is  obtained 
from  (6.2)  ; 

i     z  , 

P{Z  >  z}  =  exp  -  /   h(x)dx 
(   0         ) 


=  exp|-  prH*  Bx  +  A)dx 


=  exp 


A  £n(l  +  |)  +  |  z2  +  Az  1 


I 
=  (_£-)   exp(-  |  *-J  e 


a   ,A     ,   B   2,   -Az 


=  F1(x)  F2(z)  F3(z)  ,  (6.6) 
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where 


Vz)  =  (rfr)A 


F2(z)  =  exp(-  |  z2)  (6.7) 


F3(z)  =  e 


-Az 


All  of  the  above  are  recognized  as  being  the  complements  of 
distribution  functions.   In  effect,  the  distribution  of   Z 
is  that  of  the  minimum  of  three  independent  random  variables: 


P{Z  >  z}  =  P{X,  >  z}»P{X2  >  z}'P{X   >  z}  ,       (6.8) 


X.   having  the  distribution   F.  =  1  -  F.  (i  =  1,2,3).   This  fact 
leads  directly  to  an  easy  procedure  for  simulation  of   Z   by 
simply  obtaining  the  smallest  from  among  the  realization  of 
X,  ,  X2 ,  and   X-..   The  advantage  of  the  previous  method,  based  on 
(3.6)  for  instance,  is  that  only  one  realization — that  of  an 
exponential  in  that  specific  example — leads  to  the  realization 
of   Z.   This  is  not  only  computationally  attractive,  but  seems 
to  facilitate  the  application  of  such  Monte  Carlo  variance 
reduction  techniques  as  control  and  antithetic  variables 
cf .  Hammersley  and  Handscomb  (1964) . 
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VII.   OBTAINING  SPECIFIED  HAZARD  BEHAVIOR  BY  SIMPLE  SAMPLING 

The  development  of  the  last  section  illustrates  one  manner 
in  which  hazard  behavior  may  be  conveniently  represented  and 
simulated.   We  now  show  how  such  behavior  may  alternatively  be 
obtained  by  simple  simulation,  i.e.  from  one  realization  of  a 
basic  (possibly  exponential)  random  variable. 
Refer  to  (3.5) ,  in  which 

Z  =  G(X)  (7.1) 

and,  if   G(-)   is  monotonically  increasing, 

z(p)  =  G(x(p))  ,  (7.2) 

z (p)   and   x(p)   being  the  p*100%  percentiles  of   Z   and   X 
respectively.   Then  the  counterpart  to  (4.14)  that  results  from 
differentiation  of  (7.2)  is  the  expression 


hz(z(P))  =  hx(x(P))  G,u\p))  =  hx(x(P))  jz^r      (7-3) 


Consequently,  if  one  specifies   h  (z)   as  a  suitable  function  of 

the  "time"   z,  and  specifies  the  distribution  of  the  stochastic 

variable   X — and  hence  its  hazard,  h  — there  results  a  differential 

x 

equation  for   z (x)  =  G(x): 

h  (z)  ~  =  h  (x)  ;  (7.4) 

z     dx     x 
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integration  then  provides  the  desired  transformation,  G.   In 
other  words,  we  seek   z (x)   satisfying 


z  x 

/   h  (u)du  =  /   h  (v)dv,  (7.5) 

0    z        0    x 


which  can  sometimes  be  carried  out  in  a  useful  closed  form. 


Example  7.1.   Refer  to  the  example  of  Section  VI,  wherein  h 
is  given  by  the  expression  (6.2)  and  we  assume  that   X   is 
exponential,  so   h    is  constant.   Then  in  order  to  determine 


G(x)  =  z (x) ,  solve  the  equation 


/    — § —  +  Bv  +  X 1  dv  =  /   du 
n   Lv  +  a         J      J 


or 


A  ln(l   +   -)  +  |  z2  +  Az  =  x  (7.6) 

a     2 


Closed-form  solution  of  this  expression  for   z   in  terms  of   x 
is  of  course  impossible.   One  possible  approach  is  purely 
numerical:   find  an  approximate  solution,  zfl (x) ,  e.g.  the  appro- 
priate solution  of  the  quadratic 

|  z2  +  Xz  =  x  (7.7) 
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and  then  correct  the  result  by  a  few  Newton-Raphson  iterations 
In  other  words,  put 


w  .  -  A  +  A2  +  2BXX_  .  (78) 


now  apply  Newton  to  obtain  an  improved  solution 


A  ln(l    +    z  /a) 

Z~  (x)   =  Z,   -  — n ; 1 ; = ; r-  (7.9) 

2        1    A/(z,  +  a)  +  Bz,  +  A 


which  will  be  feasible  if   0  <  z_  .   The  process  can  be  iterated 
(the  numerator  will  change  after  the  first  iteration) .   If  one 
wishes  to  use  this  model  it  may  actually  be  desirable  to  start  by 
solving 

2-  z2  +  (A  +  -)z  -  x  =  0  (7.10) 

2  a 


for   z. ,  in  which  case  the  numerator  will  not  be  as  shown  in  (7.9) ; 
convergence  may  be  more  rapid. 


Example  7.2.   Change  the  hazard  representation  of  the  previous 
example  as  follows:   let 


h  (v)  = 7   +  Bv  +  A;   (A,  a,  B,  A  >  0)  (7.11) 

Z       (v  +  a) 


then 


/  h  (v)dv  =  —nr^ — r  +  |  z2  +  Az  (7.12) 

^    z   ;      a (z  +  a)    2 
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Now  the  solution  of  the  equation 


+  ?  z2  +  Az  =  x  (7.13) 


a (z  +  a)    2 

can  be  carried  out  in  closed  (if  cumbersome)  form,  since  (7.13) 

is  a  cubic.   Hence  an  explicit  representation  of   z (p)   in  terms 

of  an  exponential's   x(p)  can  be  provided.   Perhaps  more  importantly, 

a  direct,  simple,  simulation  of  a bathtubbed  random  variable  is  at 

hand. 


VIII.   CONCLUSION 

The  purpose  of  this  report  has  been  to  show  that  rather 
complex  distributional  behavior  can  be  represented  in  terms  of 
simple  standard  distributions.   We  have  concentrated  here  upon 
alterations  of  the  exponential,  but  it  is  obvious  that  other 
distributions  can  be  treated  similarly. 

In  a  later  report  we  will  illustrate  the  use  of  this  technique 
in  data  analysis,  and  in  the  simulation  of  non-homogeneous 
stochastic  point  processes. 
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